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ABSTRACT 

A model of the gravitationally evolved dark matter distribution, in the Eulerian space, 
is developed. It is a simple extension of the excursion set model that is commonly used 
to estimate the mass function of collapsed dark matter haloes. In addition to describing 
the evolution of the dark matter itself, the model allows one to describe the evolution 
of the Eulerian space distribution of the haloes. It can also be used to describe density 
profiles, on scales larger than the virial radius, of these haloes, and to quantify the way 
in which matter flows in and out of Eulerian cells. When the initial Lagrangian space 
distribution is white noise Gaussian, the model suggests that the Inverse Gaussian 
distribution should provide a reasonably good approximation to the evolved Eulerian 
density field, in agreement with numerical simulations. Application of this model to 
clustering from more general Gaussian initial conditions is discussed at the end. 



Key v^rords: methods: analytical - galaxies: clusters: general 
cosmology: theory - dark matter. 
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1 INTRODUCTION 

The hypothesis that, in comoving coordinates, initially 
denser regions contract more rapidly than less dense regions, 

■ and that sufficiently underdense regions expand, is simple, 

■ reasonable, and powerful. As a consequence of this expan- 
, sion and contraction, the density distribution in the initial 
' Lagrangian space will be different from that in the evolved, 
' Eulerian space. Suppose that, as the universe evolves, the 

number of expanding and contracting regions is conserved — 
only their comoving size changes — and the mass within each 
such region is also conserved. If we have a model for the way 
in which the evolution of the size of a region depends on its 
initial size and density, and we also have a model for the 
initial number of regions as a function of initial size and 
density, then we can compute the distribution of sizes and 
densities at some later time. For example, suppose that the 
initial Lagrangian density distribution is a Gaussian random 
field, and that the evolution of regions is given by the spher- 
ical collapse model. Then it should be possible to construct 
a model for p{Alo\R, z), where p{Mo\R, z) is the fraction of 
regions of size R that, at z, contain mass Mq. The quantity 
p{Mo\R,z) is often called the Eulerian probability distribu- 
tion function. 

In the Press-Schechter (1974) approach, at any given 
time, all matter is in the form of collapsed objects, usually 
called haloes, and the distribution of halo masses evolves 
with time. At any time, the matter within a randomly 



placed cell 7? is divided among many collapsed haloes. Thus, 
p{Mo\R, z) depends both on the halo mass function, and on 
the spatial distribution of the haloes. Bond et al. (1991) 
showed how to estimate the evolution of the halo mass func- 
tion if the initial Lagrangian distribution is Gaussian (also 
see Lacey & Cole 1993) . They did not show how to estimate 
the spatial distribution of these haloes, but Mo & White 
(1996) showed how this might be accomplished. This paper 
develops a model that combines, self-consistently, the Bond 
et al. (1991) excursion set approach with the Mo & White 
(1996) model for the Eulerian space halo distribution. The 
model developed here allows one to simultaneously describe 
both the distribution of the dark matter, i.e. p{Mo\R,z), 
and that of the haloes. See, e.g.. Mo & White (1996) for 
why such a model is useful. 

Section |^ describes the model. It shows why the dis- 
tribution of first crossings of a barrier whose height is 
not constant, by Brownian motion random walks, is use- 
ful. The shape of the barrier associa ted with the spheri- 
cal collapse model is given in Section |2.2| . The relation be- 
tween the first crossing distribution and the Eulerian dis- 



2.4 



tribution p[Mo\R,z) is discussed in Section Section 
shows that, in the context of the model developed here, the 
halo mass function is related to the small cell size limit (i.e. 
R ^ 0) oi p{Mo\R,z). It discusses the Bond et al. (1991) 
excursion set results in this context, and then shows how 
the model can be used to describe the spatial distribution 
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of the haloes as well. Th e results of M o & White (1996) are 
Section 



discussed in Section 2.5 



2.6 shows that the associ- 



ated two barrier problem can be used to provide information 
about the evolved density profile, and also about the way in 
which matter flows in and out of Eulerian cells. 

The first crossing distribution associated with the 
spherical collapse barrier must be obtained numerically. 
Therefore, to illustrate the usefulness of our approach. Sec- 
tion 1^ shows the results of assuming that the initial distri- 
bution is white-noise Gaussian, and that the barrier shape 
is simpler than that associated with the spherical collapse 
model. In Section ^, the barrier is assumed to be linear for 
a number of reasons. Firstly, this linear barrier can be un- 
derstood as arising from a simple variant of the spherical 
collapse model (Section [3.l| ). Secondly, the barrier crossing 
distribution can be computed analytically (e.g. Schrodinger 
1915). The details of the derivation are presented in Ap- 
pendix Thirdly, the Eulerian probability distribution as- 
sociated with the first crossings of this linear barrier is In- 
verse Gaussian (Section 3.3), and the Inverse Gaussian pro- 
vides a good fit to the Eulerian distribution measured in 
numerical simulatio ns o f clustering from white noise initial 
conditions (Section 3.6). Finally, studies of clustering from 
Poisson initial conditions also suggest that this barrier shape 
is a useful approximation (Sheth 1998). 

The Bond et al. (1991) results are derived, within the 
context of this linear barrier model, in Section Ana- 
lytic expressions for the Eulerian distribution function as 
well as the halo-mass and halo-halo correlations are de- 



and 3.4 



rived for all scales, and for all times, in Sections 3.3 
Section provides various formulae associated with the 
two-linear-barrier problem. A final section summarizes the 
results. It also discusses how the model can be extended to 
describe clustering from more general initial conditions. 



2 THE EXCURSION SET MODEL 

The first subsection defines various Eulerian space quan- 
tities of interest in this paper. Since these definitions are 
standard (e.g. Peebles 1980), no further references are given. 
Subsequent subsections develop a model which allows one to 
estimate the evolution of these quantities. 



2.1 The matter and halo distribution in Eulerian 
space 

Imagine partitioning the Eulerian space Kot into a large 
number of cells, each of size V = (47r/3)-R^, at time z. Since 
the volume of each cell is V, the total number of such cells 
is A'tot = Vtot/V. Let p{Mo\R, z) AMq denote the fraction of 
these cells that contain mass between Mo and Mo -I- dMo. 
Then, at a given z, 

p(Mo|i^,.)dMo= ^(^y = ^^(^;^i^)dMo, (1) 



where we have not bothered to write z explicitly on the 
right hand side. Define A = (1 + 5) = Mo/pV , where p is 
the average comoving density in the Eulerian space. That 



IS, 



TT^^T^ MoN{Mo\R,z)dMo. 

Vtot Vtot J Q 



(2) 



Then 

p(A|i?, z) dA = p{Mo\R, z) dMo, 



(3) 



and p(Ai_R, z) is the probability distribution function of the 
density in Eulerian space at z. Since dA = AMo/{pV), 



p{A\R,z)dA^ Ap{A\R,z) dA=l 



(4) 



for all z. 

Let n(A/i, Sd), where Sd = (5co(l + zi) and Sco is some 
constant that will be determined later, denote the average 
number density of Mi haloes identified at zi. On average, 
the number of such haloes within an Eulerian cell V is 
n{Mi, Sci)V . Let A''i denote the number of such Mi haloes 
within the ith Eulerian cell. Suppose we classify all Eule- 
rian cells by the mass Mo within them at some z < zi. Let 
A'^(1|0) denote {Ni}q, where the average is over only those 
Eulerian cells that contain mass Mo at z. Clearly, 



(5) 



(6) 



n{Mi,5d)V^ / N{1\0) p{A\R,z)dA, 
Jo 

since equation (^) shows that p{Mo) dAdo = p(A) dA. 
Define 



5h(l|0) = - 



^(1|0) 
n(Mi,4i)V 



This is the number, relative to the average number of 
(Ml, 6ci) haloes in Eulerian cells V, of such haloes that are 
within Eulerian cells which contain mass Mo at z, minus 
one. The cross correlation ^hm(Mi, (5ci |-R, z) between haloes 
and mass in Eulerian space, averaged over Eulerian cells of 
comoving size V, at the epoch z, is 



^bU^\R,z) 



4(110) 5 



\A-i)^mm3AiA, 

^ ' n{M,,Sd)V 



(7) 



where we have used equation (^ to set (A) = (1 + 5) = 1, 
and so (S) — 0. 

Let ^hh(Mi, M2, 5ci\R, z) denote the (volume average of 
the) correlation function of Mi and M2 haloes identified at 
the epoch Zi, averaged over all comoving Eulerian cells V, 
at the epoch z < z\. This average can be computed in two 
steps. Let A''i and N2 denote the number of Mi and M2 
haloes within an Eulerian cell V . Classify all Eulerian cells 
by the mass Mo within them at z. Let C(12|0) denote the 
average over all Mo cells of {N1N2): C{12\G) = {NiN2)o- 
Then 



l + 6.h(12|i?,z) 



C(12|0) piA\R,z) 
'^n{Mi,6ci)VniM2,5ci)V 



dA. 



(8) 



Higher order moments of the halo distribution can be defined 
similarly. The remainder of this section develops a barrier- 
crossing, excursion set model which allows one to estimate 
all these Eulerian quantities. 



2.2 The spherical collapse model 

We will assume that the total mass and comoving volume 
of the evolved Eulerian space is the same as that initially 
(i.e., in the Lagrangian space). Then the average density in 
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the two spaces is the same: po ~ p, where here and below 
quantities with subscript zero are in the Lagrangian space. 

We wiU also assume that the initial overdensity fluctua- 
tions in the Lagrangian space are small: So <^ 1 initially. If p 
is the average background density, then Mo = poVb(l + (So) ~ 
pVo, where Vb = (47r/3) Rq. That is, the initial mass Mo, vol- 
ume Vo and size Ro are all equivalent variables. Consider a 
region that initially contains mass Mq and has initial over- 
density So- At some later time, it has size R{z), so that the 
overdensity within it is 5, where 

^3 



1 + 5: 



Mo/pV = {Ro/RY 



(9) 



Notice that A is the same Eulerian quantity of the previous 
subsection, and that it is simply the ratio of the Lagrangian 
volume to the Eulerian volume of the region. 

The spherical collapse model (e.g. Peebles 1980) allows 
one to describe the evolution of such a region. In particular, 
it provides relations between the initial size Ro, the initial 
overdensity 5o, the time z, and the evolved Eulerian size at 
that time R{z). If -Ro, ^o, and z are given, then R{Ro,5q, z) 
is determined by the model. On the other hand, if only R 
and z are given, then the model describes a curve in the 
{Ro,So) plane: 5a{Ro\R,z). To a good approximation, this 
relation is 



5o{Ro\R,z) 
l + z 



= 1.68647 - 



1.35 

A2/3 



1.12431 ^ 0.78785 ^^^^ 



Ai/2 



/\0. 58661 



(Mo & White 1996). A simpler approximation to this rela- 
tion is 



5o{Ro\R,z) 
l + z 



5,0 A 



(11) 



(Bernardeau 1994). Whereas Bernardeau used 5co = 1.5, the 
value 1.68647 is also acceptable. 

Notice that when R ^ 0, then A — > oo, so formulae ( p^ 
and (0) both become 5o 5co (1 + z) with 5co ~ 1.686. In 
this limit 5o{Ro\R = 0, z) is independent of Ro- When R > 0, 
then 5o{Ri)\R,z) decreases monotonically as Rq decreases. 
Essentially, equation ( [ic| ) shows that a given pair R and z 
could initially have come from a range of Ro and 5o. This 
is sensible; in the spherical model, initially denser regions 
collapse more rapidly than less dense regions. Therefore, a 
region of size R at z may initially have been a small region 
containing a small overdensity, or it may have been larger 
initially, but with a correspondingly larger overdensity. 

Our model for estimating the Eulerian probability dis- 
tribution function p{Mo\R, z) works as follows. At any given 
z, imagine partitioning space into a large number of cells 
each of size V. Assume that each cell evolved according to 
the spherical model independently of the others. This means 
that the mass within each cell remains the same — only its 
comoving size changes — and the total number of cells is con- 
served. These are strong simplifying assumptions that do 
not have a rigorous physical justification. However, they al- 
low one to estimate a number of useful quantities. As with 
the Bond et al. (1991) excursion set approach, the extent 
to which the model here is able to reproduce the results of 
numerical simulations is the only real justification for these 
assumptions. 

Since each cell evolved according to the spherical model, 
and the number of such cells is conserved, to compute 
p{Mo\R, z) we simply need to specify the relative numbers of 
regions initially with (i?o, <5o) that have now evolved into re- 



gions {R, z). Clearly, this distribution depends on the initial 
distribution of fluctuations and on the evolution described 
by the spherical model. The next subsection shows how to 
do this. 



2.3 The first crossing distribution 

Consider a density field 5o{r;Qi). Recall that the subscript 
zero denotes the fact that 5o is a quantity measured in the 
Lagrangian space. Smoothing this field with a window W of 
scale Ro produces a smoothed field 5a{r;R{j): 

5o{r;Ro) = j W{\r - r'\,Ro) 5o(r' -Q) Ar' 



W{k,Ro) <5o(fc) exp(ife.r) dfc, 



(12) 



where W is the Fourier transform of the window W. By 
definition of the average density, 5o ^ as i?o — + oo, for all 
positions r. Thus, for each position in the Lagrangian space 
there is a curve, 5o{Ro), which describes the overdensity 5o 
in a window of Lagrangian size Ro centred on that position. 
Call such a curve a trajectory. The volume associated with a 
window of size Ro is Vo = Vw^^o, where Vw is some constant 
that may depend on the shape of the window, but does not 
depend on Ro- For example, Vw ~ 47r/3 for a Top Hat 
window. The mass within this volume is Mo — poVb(l -I- 5o), 
where po is a constant that denotes the average mass density 
in the Lagrangian space. Provided 5o is small. Mo ~ poVo to 
lowest order. In this approximation, there is a deterministic 
relation between the volume Vo and the mass within it: AIo = 
poVo in the Lagrangian space. 

There is also a deterministic relation between the scale 
Ro, and the variance 5*0 associated with that scale: 



5-0 = {5Ur;Ro)) = / Pik)W\k;Ro) d\, 



(13) 



where P{k) is the power spectrum of the unsmoothed La- 
grangian field. Thus, for a given P{k), the quantities 5*0, -Ro, 
Vo and Mo are equivalent. For most power spectra of inter- 
est, 5*0 increases as -Ro decreases. For what follows, assume 
that So —> oo as Ro —> 0. 

The deterministic relation between -Ro and So means 
that to each trajectory 5o{Ro) associated with a position in 
the Lagrangian space, there is a corresponding trajectory 
50(5*0). Since 5o{Ro) — > as -Ro 00, all these trajecto- 
ries start from the origin (5*0, 5o) = (0, 0). Since 5*0 —> 00 as 
-Ro — > 0, and since So is a measure of the mean square dis- 
tance of a trajectory at So from the So axis, any trajectory 
may be an arbitrary distance above or below the So axis as 
Ro 0. 

The deterministic relation between Ro and So also 
means that equation (|l^) represents a curve in the (So,(5o) 
space. Given -R and z, draw this curve, and call it the bar- 
rier -B(So|-R, z). The thick solid curve in Fig. |^ shows this 
barrier for a representative value of -R and 2. It decreases 
monotonically as So increases. For a given 2, equation (|lo|) 
shows that curves for different -R all have the same value 
B = (5co(l + z) at So = 0. Equation (0) also shows that 
-8 = when R — Ro, so curves for larger values of -R cross 
the So axis at smaller values of So. In fact, for a given z, 
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Figure 1. Examples of trajectories (thin jagged curves) traced 
out by the Lagrangian overdensity, (5q , as a function of linear vari- 
ance, Sq. The trajectories are absorbed at the barrier (thick solid 
line). Here, the barrier shape is given by the spherical collapse 
model (equation |lO[), and So oc l/Vb as it is for white noise. 



the curves for different R form a nested family (see Fig. 
This will be important below. 

Since B{So\R, z) decreases monotonically with increas- 
ing So, each Lagrangian trajectory will intersect the bar- 
rier at least once (see Fig. uh. Define the first crossing as 
the smallest value of So at which the Lagrangian trajectory 
SoiSo) intersects B{So\R, z). That is, it is the smallest So 
at which Sq = B. Bond et al. (1991) discuss why the first 
crossing distribution is so important (it solves the so-called 
cloud- in-cloud problem) . Although they only considered the 
special case in which B{So\R, z) was independent of So, their 
argument holds for the more general barrier considered here 
also. The argument is as follows. 

Since So and Mo are equivalent variables, the first cross- 
ing value of So represents the largest value of Mo at which 
the trajectory had (5o(So) = B. Since Mo and Vb are equiv- 
alent variables, the trajectory can be thought of as repre- 
senting a volume element of a Lagrangian volume Vb that 
has Eulerian size R at z. Suppose we agree that the barrier 
crossings are significant, and the mass we associate with a 
trajectory could be that associated with any of the crossings. 
If these possible assignments correspond to Mo > Mi > ■ • ■, 
then the trajectory represents concentric Lagrangian regions 
Vb > Vi > In effect. Bond et al. (1991) assume that 
successive shells never cross; whereas the actual sizes of 
Vb > Vi > • • • may change, the order of the shells is pre- 
served. This means that if the largest region Vb has Eulerian 
size V , then the subregions originally within it will remain 
within it. So the largest Lagrangian region associated with 
the trajectory, i.e., the one associated with the first crossing, 
is the one that is relevant. Counting the largest region Vb 
at once includes all the smaller ones within it, and is the 
natural way to avoid double-counting the smaller regions. 

Although Bond et al. (1991) only considered the case 
in which Lagrangian regions had collapsed to Eulerian size 
R = Q, their argument can also be applied in the R> Q case 
studied here. Recall that, for a given z, curves for differ- 
ent R are nested. Consider a trajectory 6o{So) that crosses 
the barrier B{R) for the first time at Vb. The Lagrangian 



scale Vn-i, corresponding to the nth crossing of B{R) by 
this same trajectory may be smaller than the Lagrangian 
scale Vq corresponding to the first crossing, by this trajec- 
tory, of the barrier B{R') for some Eulerian R' < R. In 
our model, such a trajectory represents a volume element of 
a Lagrangian region Vb that at z, has Eulerian size R, at 
which time it contains a subregion that originally had size 
Vo < Vb and now has size R' < R. The assumption is that 
if Vq has Eulerian size R' , and if Vn-i < Vq < Vb originally, 
then all the mass Mn-i contained within Vn^i must still be 
contained within the Eulerian region R' now occupied by 
Vq. Thus, the model developed here explicitly preserves the 
ordering that was assumed by Bond et al. (1991), and so the 
first crossing of the barrier B{So\R, z) remains the one that 
is significant. 

Let f{So,So\R,z)dSo denote the fraction of trajecto- 
ries that have their first crossing of the barrier between 
So and So + dSo. Then, for a given barrier B{So\R, z), 
f{So,So\R,z)dSo can be equated with the fraction of La- 
grangian space that is associated with regions containing 
mass Mo that each occupy an Eulerian region i? at z. If 
A'^(Mo|-R, z) dMo denotes the number of such regions, then 
each such region occupied a Lagrangian volume Vo ~ Mo/po, 
so 



f{So,5o\R,z)dSo 



VoN{Mo\R,z)dMo 

Vtot 



(14) 



where Kot = J Vo N{Mo\R, z) dMo. Thus, the distribution 
of first crossings gives an estimate of the number density 
of Lagrangian regions which contained mass Mo and had 
initial overdensity 5o, so that at z, while they still contain 
mass Mq, they occupy the Eulerian volume V . 

If we further suppose that the number of such regions is 
the same in both the Eulerian and Lagrangian spaces, then 
comparison with equation (|l]) shows that 

p{Mo\R,z)dMo = ^ f{So,5o\R,z)dSo. (15) 
Vo 

Since (Vo/V) = Mo/(pV) = (1 + (5) = A, this means that 



/(So, So\R, z) dSo = Ap(A|i?, z) dA, 



(16) 



provided that 5o is the function of So, R and z that is given 
by equation ([ic|). This shows how the barrier crossing dis- 
tribution is related to p(Aj_R, z). 

Notice that the left hand side of equation ( p^ is de- 
termined by Lagrangian space quantities (the trajectories 
which cross the barrier are Lagrangian), whereas the right 
hand side is associated with Eulerian space. A similar sort of 
relation between Lagrangian and Eulerian space quantities 
was used by Bernardeau (1994) (see the discussion at the 
start of his Section 3.2.2). The difference between his work 
and ours is that he used a Lagrangian distribution that was 
derived from perturbation theory, rather than from a barrier 
crossing model, for the left hand side of equation (p^. 

In summary: The Eulerian probability distribution 
function p{Mo\R, z) is the fraction of regions of size R that, 
at z, contain mass Mq- To estimate it, first assume that the 
initial fiuctuation So is small, so Mo/po = Vb = (47r/3)i?o. 
Next, choose a random position in the initial field. Compute 
5o in spheres 7?o centred on this position. Call the curve 
5o{Ro) centred on this position a trajectory. Given R and 
z, draw the barrier B{Ro\R, z) associated with the spherical 
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collapse model (equation |lO[). The trajectory So(Ro) may in- 
tersect the barrier B{Ro\R, z) many times. Find the largest 
value of Ro at which the trajectory intersects B{Ro\R,z). 
Call this the first crossing of B{Ro\R,z). Associate mass 
Mo with this trajectory. Since mass Mo and initial volume 
Vb are equivalent, this trajectory represents a volume ele- 
ment of a region containing mass Mo. Initially, Mo had size 
(47r/3)-Ro s-nd overdensity 5q{Ro\R,z). At z, it has size R. 
Repeat for an ensemble of such trajectories. So compute the 
distribution of first crossings of B{Ro\R, z). The fraction of 
trajectories which first cross the barrier B[Ro\R, z) at i?o 
represents the fraction of mass that is in regions of mass Mo 
that, at z, have size R. From this, the average number den- 
sity A'^(Mo|-R, z) /Vtot of such regions can be computed easily. 
Assume that the number of such regions is conserved. Then 
N{Ma\R, z)V/Vtot equals p{Mo\R, z). 

2.4 The halo distribution 

The previous subsection showed how to compute the Eu- 
lerian space distribution function using an excursion set 
model. However, the excursion set model can also be used to 
provide information about the halo distribution. This sub- 
section shows why. 

Following Bond et al. (1991), define a halo as a La- 
grangian region that has collapsed to a vanishingly small 
Eulerian size: -R = 0. Suppose z is given, and consider 
the Eulerian scale i? = 0. In the spherical collapse model, 
-B(_Ro|0, z) is independent of Ro but depends on z (equa- 
tion nCJ). Define B{R(,\0,z) = 5c{z). As before, equate the 
fraction of trajectories f{Ro,Sc) which first cross Sc{z) at Ro 
with the fraction of mass in regions Mo that, at z, have size 
R^O. Following (0), the average number density of such 
regions is n(Mo,5c) = N{Mo\0, z) /Vtot = {p/Mo) f[Ro,S,), 
and n[Mo\z) is said to be the mass function of collapsed 
haloes. This is exactly the excursion set model for the uncon- 
ditional halo mass function developed by Bond et al. (1991). 
In the present context, the unconditional mass function is 
just the Eulerian distribution function in the limit R = Q. 

Since the R = Q limit has this special interpretation, in 
what follows, it is convenient to make a distinction between 
the barrier crossing distribution in this limit, and that when 
R > Q. Below the subscript 'c' (for constant height), as 
in /c(So,(5o), denotes the barrier crossing distribution when 
R = Q, and the subscript 'R' denotes the case R> Q. 

Suppose we consider two barriers B\ = B[So\Ri, z\) 
and -Bo = B{So\R, zo), with z-i > zo- (It may help to look at 
Fig. ^.) For now, assume that R\ < R. Let /r(S'i, ^ilSo, 5o) 
denote the fraction of trajectories which first cross Bi at Si 
given that they first crossed Bo at So. Again, consider the 
limit Ri = R — for both barriers. Then Bi — Sc{zi) — Sd, 
and -Bo = ^co for all So, so both barriers have a con- 
stant height, and Bi > -Bo. Bond et al. (1991) interpret 
/c(Si, (5i|So, 5o) as representing the fraction of mass in an 
Mo-halo that was earlier in Mi-haloes. They define a con- 
ditional mass function as the average number of Mi-haloes 
that are within an A/o-halo: 7V(1|0) = (Mo/Mi) /c(l|0) if 
Ml < Mo, and A/'(1|0) = otherwise. In the context of 
this paper, /c(l|0) represents the fraction of the Lagrangian 
region Vb that has Eulerian size 7? = at «o, which was pre- 
viously in Lagrangian regions Vi that, at zi, also occupied 
vanishingly small Eulerian volumes. The reason for wording 



things in this way is that it shows how to compute other 
properties of the halo distribution in Eulerian space. 

For example, the average number of A/i haloes that col- 
lapsed at zi and are in Eulerian regions R > a,t z < zi 
can be computed as follows. Recall that (Mi, zi)-haloes are 
associated with trajectories which cross the constant bar- 
rier -Bi — Sci at Si. Since Sd > 5c(z), and since -Bo de- 
creases monotonically as So increases, each of these trajec- 
tories must have crossed -Bo — B{So\R, z) at some So < Si, 
so as to reach Sc(zi) at Si. Therefore 
nSt 

/c(Si,4i)= / /c(Si,5ci|So,<5o)/r(So,<5o) dSo, (17) 
Jo 

where <5o = B{So\R, z). Since /c(l) = (Mi/p) n(l), /c(l|0) = 
(Mi/Mo)7V^l|0), /r(0) dSo = Ap(A) dA, and A = Mo/pV, 
equation (hj) implies that 

/•oc 

n{Kh,5d)V= / N[Kh,5d\Mo,5o)p{Mo\R,z)AMo. (18) 
Jo 

If we set 

Ar(l|0) =AA(Mi,4i|Mo,5o), (19) 

where -/V(1|0) was defined in Section ^, just prior to equa- 
tion (|^), and &o in A/'(1|0) is given by equation (^^, and we 
recall that p(Mo)dMo = p(A) dA jfequation |^^then this 
expression is the same as equation (g). Equation (|l^) shows 
that the average number of Mi-haloes in Eulerian cells V 
which contain mass Mo at z is equal to that in those La- 
grangian regions Vb which had initial density 5o given by 
equation (px|). 

The main reason for writing this out explicitly is that it 
shows how statistics in the Lagrangian space can be used to 
compute statistics in Eulerian cells V . Namely, the assumed 
conservation of number of regions in the two spaces (equa- 
tion |l^ equals equation ^ means that an average over Eu- 
lerian regions of size V is equivalent to averaging over those 
Lagrangian regions that have evolved into Eulerian regions 
of size V (equations ^ and [l9| ). In practice, this means that 
Eulerian space quantities like ^hm and ^hh can be computed 
simply by using the appropriate value for So (that given by 
equation |l^ in the associated Lagrangian space quantities, 
and then summing over Mo , weighting the contribution from 
each Mo appropriately. This weighting is given by the barrier 
crossing algorithm, because the barrier crossing distribution 
/r(So,(5o), with So given by equation (p^, represents the 
fraction of Lagrangian space that is associated with regions 
of mass Mo that, because they had initial overdensity So, 
are Eulerian regions of size V at z. 

Section ^ illustrates how this works. It presents a model 
in which all these quantities can be computed analytically. 

2.5 Relation to the work of Mo & White 

Before moving on, it is interesting to compare this model 
with previous work. Our equation ( p^ is exactly the same as 
that assumed by Mo & White (1996). They assumed that the 
Eulerian space quantity -/V(1|0) could be computed by sim- 
ply substituting the spherical collapse value of So{Ro\R,z) 
into the Lagrangian formula A/'(1|0). They showed that the 
resulting formula for -/V(1|0) provided good fits to the cor- 
responding Eulerian space quantity measured in their nu- 
merical simulations of clustering. Since, in this regard, our 
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model leads to the same formula as that of Mo & White, the 
agreement with simulations provides some justification for 
the strong simplifying assumptions which lead to our model. 

Our model extends their results in the following way. 
Our equation ( p^ ) shows that, to compute averages in Eu- 
lerian space, knowledge of the Eulerian space distribution 
function p{A\R,z) is necessary to then take the correct av- 
erage over Lagrangian regions Vo that have become Eule- 
rian regions V. Mo & White also knew this, but they did 
not know how to compute p{A\R,z). Therefore, they as- 
sumed that it could be taken directly from their numerical 
simulations. Ifowever, as our approach shows, p{A\R, z) is 
related to the shape of the boundary (cf. our equation nffl . 
This relation must not be ignored. To see why, suppose (as 
Mo & White did) that a Lognormal can be substituted for 
p{A\R, z). If one uses the spherical model for 5o{Ro\R, z) in 
A/"(liO) and then does the integral on the right hand side 
of (lla), one finds that it does not equal the correct value 
for the left hand side. This is because, as our approach 
shows, self-consistency requires that p{A\R,z) depend on 
the boundary shape; once the boundary is specified, one 
is no longer free to choose any arbitrary distribution for 
p{A\R,z). This may partially explain why the Mo & White 
results using a Lognormal are somewhat worse than when 
they use the distribution measured directly in their simula- 
tions. 

Since our approach shows how to derive p{A\R, z) from 
the boundary shape, our approach can be thought of as de- 
riving, self-consistently, the Mo & White (1996) spherical 
evolution formulae by a simple extension of the Bond et al. 
(1991) excursion set approach. 



2.6 The two barrier problem 

The previous subsection showed that problems involving the 
crossing of two barriers, in which one barrier was assumed 
to be constant, and the other not (i.e. the Eulerian scale was 
i? = for one barrier, but R > Q for the other), could be 
used to provide information about the distribution of haloes 
within Eulerian cells V. The case in which J? > for both 
barriers is also interesting. 

Consider two barriers B{Rii\R, Z) and b{ro\r,z). If z — 
Z and r < R, then the two barrier problem is associated 
with the joint distribution of the mass within two concen- 
tric Eulerian cells. Therefore, it may be useful for estimating 
the evolved Eulerian space density profile. Of course, in the 
model, haloes collapse to zero radius, and this is unreal- 
istic. In practice, haloes virialize at some fraction of their 
turnaround radius. Therefore, the solution of this two bar- 
rier problem is probably only useful on scales larger than 
the virial radius. 

If, instead, r — R, but z > Z, then the two barrier 
problem describes the matter within the same comoving cell 
at two different epochs. That is, it can be used to quantify 
the way in which matter flows in and out of (comoving) 
Eulerian cells. 




t/tjsc) 

Figure 2. The physical radius of a perturbation in units of the 
spherical model turnaround radius, as a function of time in units 
of the spherical collapse model turnaround time. The solid curves 
show the spherical model, and dashed curves show the linear bar- 
rier model studied in this paper. The two curves for each line type 
are for denser perturbations (which recoUapse) and under-dense 
perturbations (which do not). 



3 THE LINEAR BARRIER AND 

WHITE-NOISE INITIAL CONDITIONS 

This section shows how the model described in Section]^ can 
be used. It assumes that the initial Lagrangian space distri- 
bution is white-noise Gaussian and studies the first crossing 
distribution of a barrier whose height decreases linearly as 
Vo — (47r/3)-Ro decreases: 



B{So\R,z) =5,o{l + z) 



5co (1 + z) 



(20) 



where So = 1/poVo, and A = Mo/pV = (Ro/Rf- 

This barrier shape was chosen for a number of reasons. 
First, this shape is a simple approximation to the spherical 
collapse barrier (compare equation HOl and see Fig. 0) , and 



3.1 



the associated dynamical evolution is similar. Section 
shows how the evolution of a spherical perturbation in this 
model differs from that in the usual spherical collapse model. 
Second, the first crossing distribution of a linear barrier is 
known, so the associated Eulerian distribution can be writ- 
ten analytically; it is Inverse Gaussian. The Inverse Gaus- 
sian is a good approximation to the distribution measured in 
numerical simulatio ns o f clustering from white noise initial 
conditions (Section 3.6), so the model may also be realis- 
tic. Finally, for this barrier, the halo-mass and halo-halo 
correlations can all be computed analytically. 



3.1 Relation to the spherical collapse model 

The linear barrier studied in this section is associated with 
the following model for the collapse of objects. 

Let R{z) denote the comoving size of an object at the 
epoch z. Then R{z) = Rq initially. If the object is in an 
underdense region, then its comoving size will increase, else 
it will decrease. Trajectories with extrapolated linear over- 
density So greater than 5co are associated with collapsed 
objects. Collapsed objects have R{z) = 0. Equation ( |2o| ) 
implies that, before collapse. 
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ScO 



(21) p{So,So\S',5')dSo 



dSo 



The radius of an object in proper, physical coordinates 
is Rp{z) — R{z)/{1 + z). Objects which collapse have a 
turnaround radius — the maximum value that Rp{z) attains. 
This occurs at 



(l + Zta) = 

at which time 
R{zt.) _ 1 



Ro 



41/3 ■ 



(22) 



(23) 



In comparison, collapsing objects in the spherical model 
reach turnaround at 



(l + .ta)=4^/^f , 
OcO 

at which time 

-R(zta) _ (1 + 2ta) -Rp(«ta) 



6 4^/^ _ M \ 



2/3 



(24) 



(25) 



Ro Ro 

(e.g. Peebles 1980). Fig. ^ shows that, for overdense per- 
turbations, turnaround in this model (dashed lines) occurs 
later, and at a larger radius, than in the spherical model 
(sohd lines). To compensate, underdense regions expand less 
rapidly in this model than in the spherical model. 

3.2 The constant barrier and statistics in 
Lagrangian space 

In the limit R 0, equation (^o|) for the barrier shape 
becomes independent of ^o, as it does for the spherical col- 
lapse model. In this limit, the first crossing distribution is 
known (e.g. Bond et al. 1991), although we give a different 
derivation of it below. Therefore, to set notation, this sec- 
tion summarizes various known results about the crossing of 
a constant barrier by the Lagrangian trajectories described 
earlier, under the assumption that 5o{r;0) is a Gaussian 
random field. 

Let p{Sq,5o) dSo denote the fraction of trajectories 
that have value between So and So + dSo at So- Also, let 
p(5'o, ^olS", 5', first) d(5o denote the fraction of those trajec- 
tories which first crossed the barrier B{So\R, z) between S' 
and S' + dS' , which have value between So and So + dSo at 
So>S'. Then 
/■So 

p{So,5o)^ p{So, So\S', S', &r8t) f{S',5'\R,z) dS', (26) 



provided that at So, So > B{So\R, z). If both p{So,So) and 
p{So, So\S' , 5' , first) are known, then this is a Volterra equa- 
tion of the first kind, so it can be solved numerically, by 
forward substitution, to yield f{S',S'\R,z). 

Suppose that the probability distribution of the density 
in the Lagrangian space is Gaussian. Then 



p{So, So) dSo 



dSo 
V2^, 



exp 



ii. 

2So 



If the window W is sharp in fc-space, then 
p{So,So\S',S', first) dSo=p{So,So\S',S') dSo, 
and 



(27) 



(28) 



^27r(5o - 5") 



exp 



(^0 - S')^ 
'2{So~S') 



(29) 



(e.g. Bond et al. 1991). Thus, if the Lagrangian space dis- 
tribution is Gaussian, and the window function is sharp in 
fc-space, then the distribution of first crossing times is easy 
to compute, whatever the boundary shape. Equations ( |28[ ) 
and ( p9| ) are also correct for a Top Hat filter in real space, 
if the Gaussian field is white noise. So, the distribution of 
first crossing times is easy to compute for white noise also. 
In either of these cases, a trajectory So{So) resembles the 
motion of a particle undergoing standard Brownian motion 
with zero drift. Bond et al. (1991) and Lacey & Cole (1993) 
used this fact to compute f{So,So). 

Notice, however, that when the boundary is a constant, 
B{So\R = 0,z) = (5co(l + z) = Sciz), then the form of 
/(So, So) can be obtained directly from equation (p^). First, 
take the derivative with respect to So on both sides of equa- 
tion (p^), and evaluate it at So = Sc{z). Next, set S' — Sc{z). 
Then the integrand on the right hand side is zero, except 
when S' = So- Thus, equation ( p6| ) implies that 



/c(So,5cN) = 



Scjz) 

So 

Sciz 



p{So,Sc[z]) 



^2^ 



exp 



Sjiz) 
2So 



dSo 
So ' 



(30) 



where the subscript 'c' denotes the fact that this is the dis- 
tribution of first crossings of a constant boundary. 

If the trajectory is known to start from (So, 5co), rather 
than from the origin, then the distribution of first crossings 
of the constant barrier Sci — 5co(l + zi) can be solved simi- 
larly. Provided Sd > Sco and Si > So, 

/c(Si,5ci|So,(5co) = ( - S ° ) 'P{S\,Scx\So,Sco)- (31) 

These expressions for the first crossing distribution are the 
same as those derived by Bond et al. (1991) and Lacey & 
Cole (1993). 

Following Bond et al. (1991), treat the parameter S^iz) 
as a pseudo-time variable; it decreases as the universe 
evolves. Then associate the first crossing distribution of the 
barrier S^iz) with a number density of regions containing 
mass Mo in the Lagrangian space, and call every such re- 
gion a dark matter halo. Then the number density of M\ 
haloes at the epoch labelled by S^i = Sc{z\) is 

n(Afi,4i) dMi = ^ /c(Si,4i) dSi. (32) 
All 

This is consistent with equation (p^, and n{Mi,Sci) is 
sometimes called the unconstrained halo mass function at 
the epoch zi- 

Similarly, the average number of (Mi, 5ci)-haloes that 
are within an (Mo, 5co)-halo is 



AA(1|0).(^) /.(Si,4ilSo,4o)^ 



(33) 



(e.g. Lacey & Cole 1993). This is the constrained mass func- 
tion. 

The joint distribution of the number A'^i of {Mi,Sci)- 
haloes and the number N2 of (M2, (5ci)-haloes, that are both 
within the same (Mo, 5co)-halo, averaged over all (Mo,(5co)- 
haloes is 
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C(12|0) = (iViiV2,5ci|Mo,5co) = AA(1|0) AA(2jlO), (34) 
where A/'(1|0) is given by equation (|33|), and 



Af{2\10)dM2 = fc{S2,S,i\Soi,5oi) dSa (35) 

(Sheth 1996b). Here Soi = S(Afo-Mi) = (Mo-Mi)-i and 
Soi is the overdensity in the remaining volume Vb — Vi that 
is not occupied by the A^i-halo. That is, 

Mo - Ml 



1 + So 



Po{Vo-Vi)' 



(36) 



However, Mo = poVb(l + 5co) and A/i = poVi{l + Sd), so 
that, to lowest order, 



Sci — Soi = {Sci — (5co) 
~ ((5ci — 5co) 



Mo(l + (5ci) 



Mo(l + (5ci)-Afi(l + <5co) 
Mo 



(37) 



Mo - Ml 

Most of the formulae in this section are not new, but 
they will all be useful below. They have been included to 
set notation, and because, as the previous section showed, 
if these expressions are known, then simple transformations 
allow one to compute the associated Eulerian space quanti- 
ties. 



3.3 The linear barrier-crossing distribution and 
the matter distribution in Eulerian space 

The previous subsection described the first crossing distribu- 
tion associated with a constant barrier by trajectories asso- 
ciated with a Lagrangian field that is Gaussian white noise. 
Recall that the constant barrier is got from the linear barrier 
of equation ( |20[ ) by setting R = 0. When R > 0, then direct 
substitution shows that 



/r(5'o,(5o) d^o 



B(0\R,z) 



exp 



B^{So\R,z)\ dSo 



2So 



So 



where B is given by equation (|20[), satisfies equation 
The identity 



fn{So,5o) dSo 



,(38) 



(39) 



which follows from the definition of the modified Bessel func- 
tion of the third kind, is useful in proving this result. A full 
derivation of this distribution is given in Appendix 

Equation ( p8| ) is known as the Inverse Gaussian dis- 
tribution. It was first discovered in the context of Brownian 
motion by Schrodinger (1915). Folks & Chhikara (1978) give 
a review of its role in other fields. Here, the distribution rep- 
resents the fraction of Lagrangian space which is associated 
with regions of mass Mo (So) provided the initial distribution 
is white-noise Gaussian, and the dynamics are a simple mod- 
ification of the usual spherical collapse model (as described 
in the previous subsection.) 

For white noise, A/q — 1/So, so, for white noise, equa- 
tion (^) in ( |l6| ) implies that the Eulerian distribution func- 
tion that is associated with the linear barrier is 



p(A|i?)dA = pVSofRiSo,So) dSo 



exp 



2CmA 



dA 



(40) 



where 

^n.{R)=iSc0{^ + 4'pVr 



and p — 1, 



(41) 



and we have not bothered to write the z dependence of 
p{A\R,z) explicitly. Notice that p{A\R) is also an Inverse 
Gaussian distribution. The mean of this distribution is jj, — 
1, and the variance is ^m. The higher order moments of this 
distribution are simple. If 5„ is the nth cumulant of this 
distribution, then 



(2ri-3)!! 



(42) 



and it is independent of cell size V. This is simply a conse- 
quence of the fact that, in the model, haloes have vanishingly 
small Eulerian sizes: R = 0. 

As 0, this Inverse Gaussian tends to the Gaus- 

sian distribution. For fixed V, this happens if 5co(l + z) ^ 1. 
Since z is a pseudo-time variable, this means that, at suf- 
ficiently early times, the Eulerian distribution is Gaussian 
for almost all scales V. This is sensible, since at early times, 
the Eulerian and Lagrangian distributions should be similar, 
and, by hypothesis, the Lagrangian distribution is Gaussian. 
For a given z, ^ as 1/ 3> 1. This means that, even at 
later times, the Eulerian distribution function, measured in 
sufficiently large cells V, appears Gaussian. 

When V 0, then the term in the exponent of equa- 
tion @) tends to A/^m = 5?o(l + ^f/So, so p{A\R) V 
times equation (^2|), which is the same as 1/A times equa- 
tion (js^). This is consistent with the fact that, in this limit, 
equation (^^ shows that the barrier shape tends to the con- 
stant (5co(l + z), and the first crossing distribution of a con- 
stant barrier is given by equation (po|). This shows explicitly 
that the unconstrained mass function of ( ^^ is essentially 
the same as the Eulerian distribution function, in the limit 
of vanishing cell size V. All these asymptotic relations are 
sensible. 



3.4 The halo distribution 

A little algebra shows that when the barrier shape is lin- 
ear, then equation (^^ is satisfied for all R. Recall that 
that expression simply expresses the fact that the fraction 
of trajectories which first cross the constant barrier Sd at 
5*1 is equal to the fraction of trajectories which first cross 
the linear barrier (associated with the Eulerian scale R) 
at So = i3(So|-R, 2o), and then cross the constant barrier 
Sci > Sco at Si, summed over all So < Si. Of course, this im- 
plies that equation ( |l8|) i s also satisfied: n{Mi, Sci)V is equal 
to V times equation (|32|) , as it should. This shows explicitly 
that, for the linear barrier, the number density of haloes in 
the Eulerian space is equal to that in the Lagrangian space, 
and this number density is just what is required by the Bond 
et al. (1991) excursion set approach. 

The cross correlation between haloes and mass can be 
worked out similarly. Start with equation (^). Use equa- 
tion ( p^ ) to write p{A\R, Zohs) in terms of the Lagrangian 
barrier crossing distribution /r(So,(5o)- Then use ( p^ ) for 
fR{S(uSo), and recall that So — B{So\R, Zohs) where equa- 
tion (pd) gives B. This means that B{0\R,z) = 5co(l + -Zobs), 
where Zobs < zi. Then use equation (h9|) to set Ai'(l|0) = 
AA(llO), and use (|i) for Af{l\0) with So = B{So\R, Zobs)- 
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With these substitutions, the integral in equation 
be solved analytically: 



^.^m) - ^ (^) + 

pV \l + Zobs / 



A/*(Zobs) / Zl 



pV 



where 



Sci = <Sco (1 + Zl), S^o/So = Mo/M*o, 

M.(Zobs) = A/.o (1 + Zobs)"^ 



and 



(43) 



(44) 



Notice that when zi — Zohs, then ^hm ~ Mi/pV. That is, 
the only correlation that arises is that which is due to the 
fact that there is a halo of mass Mi inside V, so there is 
certainly mass Mi inside V. Since 

1 _ M,(2obs) 



^■"^■^^ 52^(1 + 2obs)2py pV 
equations (^ and (g) imply that 

^hm{l\R) ^ Ml f 1+Zl 

Cm(i?) A/.(Zobs) 



/ 1 + Zl \ ^ / Zl - Zobs \ 
U+Zobs/ V 1 + Zl / 



(45) 



This ratio is independent of the Eulerian cell size V . Again, 
this a consequence of the fact that the model assumes that 
haloes have zero volume, so if a halo is within a cell, then 
all its associated mass is also. 

This expression has an interesting relation to previous 
work. Mo & White (1996) argue that, when \Sco\ < Sci, 
then, in the limit of large V, 



(5h(l|0) S)^ 



{vl - 1)/Jc0 

5cl/[5c0(l + Zobs)] 



(46) 



where vi = S^i/Si. Now, 5^i/5co = (1 + zi), and {5"^) ^ 
for large V , so this expression is the same as what is required 
by equations ( ^^ and ( ^ ) if 5co = 1. The extra factor of Sco 
arises from the fact that, for large V the density fluctuation 
S in most Eulerian cells is small, soA=l+<5~l for most 
V. In this limit, equation (^o|) shows that 5o/(l + z ) ~ 5 in 
the collapse spherical model, whereas equation (bd) shows 
that it is ~ 5cO & for the linear barrier. 

Let Chh(Mi , M2, 5ci|J?) denote the correlation function 
of Ml and M2 haloes identified at the epoch z\, averaged 
over all comoving Eulerian cells V, at the epoch Zobs < zi. 
Recall that this average can be computed in two steps. All 
Eulerian cells can be classified by the mass Mo within them. 
Further classify all Eulerian cells which contain mass A/q 
by their Lagrangian overdensity Sq. Since the number of 
such regions is conserved, the joint distribution of the num- 
ber of Ml and M2 haloes, N1N2, averaged over all such 
(Mo, (5o)-regions, is the same in both spaces, provided we 
set So = B{So\R, z). This average gives C(12|0) of equa- 
tion @, where So = B{So\R,z) and B is given by (|^. 
All that remains is to further average C(12|0) over the num- 
ber of such conserved regions. This is done by weighting it 
by p{A\R), integrating over A, and then dividing out the 
factors expected on average. We have gone to the trouble 
of stating things this way to show that our definition of 
5hh(12|f?) includes the effects of the scatter of halo counts 
in Eulerian regions which had the same overdensity as well 
as the same radius initially. Therefore, ^hh here denotes the 
same quantity as a^h of equation (25) in Mo & White (1996). 

Now, 1 -I- ^hh(12|-R) is given by equation (^, where 
C(12j0) is given by equation (^i]). With equations ( |l6| ) 



and ( pq ) for p(A\R), and setting So = B{So\R, z) as be- 
fore, this integral can be solved analytically. The result is 
that 



Chh(12ji^) 

^m(i?) 



Ml -I- M 
M.( 



1 +M2 /Z] 
.(Zobs) V 1 



^1 ^^obs 



/ Zl - Zobs \ ' 
V 1 + Zl J 



(47) 



where zi > Zobs, and Af, (zobs) is defined as before. Like the 
halo-mass correlation function, this ratio is independent of 
V because, in the model, haloes are point sized. Thus, in 
the model, haloes are linearly-biased tracers of the mass on 
all scales. Mo & White (1996) argue that ^hh/^m should 
equal (^hm/Cm)^ on large scales. Although our equation ( |45| ) 
for ^hm/^m is similar to theirs, the right hand side of our 
equation (j^) is not the same as the square of the right 
hand side of equation maj; they differ by a factor of [Af(l + 

Zl)/M.(2obs)(l +Zobs)P~ 

When Zl = Zobs, ^hh(12ji?) — 0, which implies that 
haloes, when they first virialize, are uncorrelated with each 
other. Further, the correlation of haloes depends on the sum 
of the halo masses, but not on the masses themselves. This 
suggests that halo-halo correlations arise because of volume 
exclusion effects only. That is, halo-halo correlations arise 
from the fact that an Mi-halo occupies a region Vi initially, 
so no other halo can occupy this region. The correlation 
function is affected by the fact that a region is included, 
but not by the exact way in which the excluded region is 
populated. 

If the halo masses Mi and M2 and their formation 
epoch Zl are fixed, then the halo-halo correlation function 
at a given comoving R increases as Zobs decreases. This in- 
crease reflects the fact that haloes which formed at zi, at 
which time they were uncorrelated with each other, must 
have merged with each other to construct the more mas- 
sive haloes present at the later epoch Zobs < zi. Suppose, 
instead, that the halo masses, the comoving scale, and the 
epoch at which the halo distribution is measured (i.e., Zobs), 
are given, and we wish to consider the effect of changing the 
halo formation epoch zi > Zobs- Equation (^^ shows that 
as Zl decreases to Zobs, Chh decreases. That is, at fixed mass, 
the halo-halo correlation function in comoving coordinates 
exhibits negative evolution. 

3.5 The two barrier problem 

The previous subsections showed how the first crossing of a 
linear barrier by Brownian walks could be used to compute 
various interesting quantities associated with gravitational 
clustering. This section shows that the associated two barrier 
problem may also be useful. In general, the linear barrier 
has two free parameters which may be thought of as the 
2/-intercept Sco and the slope /3, respectively: 



50(5-0) 



PSo. 



(48) 



Below, we consider the statistics of trajectories which first 
cross a barrier with one choice of parameters, and then cross 
another. 

First, consider two barriers 5i and ^2 which have the 
same y-intercept Sco, but different slopes, l3i and 132- Assume 
that l32 > l3i. This means that trajectories must cross the 
barrier S2 before they can cross Si . We seek an expression for 
the probability that a trajectory has its first crossing of Si at 
5*1 , given that it has its first crossing of the barrier ^2 at 5*2 < 
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Si. Recall the derivation of the barrier crossing distribution 
given in the Appendix. To a particle that started from the 
origin and has just crossed the barrier 62(32), the barrier 
Si{Si) has shape 



Sl2{Sl-S2) = SiiSl) S2{S2) 

= (132 - I3i) S2 - l3iiSi ~ S2) 



(49) 



Since this barrier is also linear, the solution is the same as 
before, except that — > 5*1 ~ ^2 and <5i2- Therefore, 
the first crossing distribution is 



fiSi,5i\S2,52)dSi 



exp 



dSi 



{I32S2 - (3iSiY 
2(Si - S2) J (Si - S2) 



(50) 



(f32-l3i)S2 
y/2Ti(Si - S2) 

Since 5cO is a pseudo-time variable, and it is the same for the 
two barriers, equation (^8|) shows that the expression above 
will be useful for computing statistics in concentric cells Vi 
and V2 > Vi, at the same epoch. Therefore, it may be useful 
for computing Eulerian density profiles within the context 
of the model. Since, in the model, haloes collapse to zero 
radius, this way of computing density profiles is probably 
only useful on scales larger than the virial radius. The virial 
radius can be derived from combining the spherical collapse 
model with the fact that the mass within the virial radius is 
given by the scale at which the trajectory first crosses Sco- 

The corresponding expression for two barriers which 
have the same slope (i, but have different values for the y- 
intercept, say, 5ci and 8^2 < 5c\ is 

(5c\ — 5c2) 



/(Si,(5i|S2,52)dSi 



X exp 



[(5c 



^27r(Si - S'2) 
5c2)-/3(5i 



5-2 



dSi 



2(Si - S2) 



{Si - &) ■ 



(51) 



Equation ( |48[ ) shows that this provides information about 
the matter in concentric cells Vi and V2 > Vi, at two dif- 
ferent epochs. Since the value of /3 is the same, the two cell 
sizes are related: (V2/V1) — {5ci/5c2)- 

More interesting, perhaps, is the case in which the cell 
sizes V are the same, but 5ci > 5c2. This case is also more 
complicated, since now it is possible to cross barrier 1 before 
crossing barrier 2. Provided 5*2 < pV , 

f{Si,5i\S2,52)ASi ^"^^ 



X exp 



^27r(Si - S'2) 

Scl2 - PljSl - S2)Y 

2{Si - S2) 



dSi 



{Si - S2) ' 



(52) 



where Sci2 ~ {Sci — i5c2)(l — /32), and /3i and /32 are given 
by using (|2^) in (ji^), with Sd and 5^2, respectively. When 
the first crossing of S2 occurs at S'2 > pV , then the trajec- 
tory will have crossed Si first. In this case, /(S2, felSi, 5i) 
is given by the expression above, provided the labels 1 and 
2 are interchanged. Bayes' rule then gives the associated 
f{Si,Si\S2,52), since 

MSi,Si) 



f{Sl,5l\S2,S2)= fiS2,S2\Sl,Sl)^ 



(53) 



MS2,52)- 

These expressions contain information about the mass in the 
same comoving Eulerian cell at different epochs. They show 



that it is possible that the mass in V first decreases and later 
increases. This shows explicitly that, in the model, matter 
can fiow in and out of Eulerian cells. 



3.6 Comparison with simulations 

This subsection shows that the Inverse Gaussian distribution 
provides a good description of the evolved distribution for 
clustering from white noise initial conditions. Before doing 
so, we first argue that such a test has, in fact, already been 
performed by others. 

The Inverse Gaussian is a limiting form of the general- 
ized Poisson distribution (GPD) 



(54) 



p(ivii^,fe) = §[7v.+iV6]"-^-^--^ 

where A'^ is the number of identical particles in a randomly 
placed cell of size R which contains iV particles on average, 
Nc = iV(l — b), and < b < 1. In the astrophysical context, 
the GPD was first derived by Saslaw & Hamilton (1984). It 
can be related to the Inverse Gaussian as follows (Fry 1985; 
Sheth 1996a). The mean of this discrete distribution is N, 
and the variance is iV/(l - bf = N{1 + N(,). Notice that as 
iV ^ 00 and fe = (1 -I- 5co)"^ ^ 1, ^ S~g, which is the 
relation required by equation (pl|). In this limit, N p{N\R, b) 
tends to the Inverse Gaussian distribution p{A\R), since 
N/N = A. 

Sheth (1998) shows that the GPD can be derived using 
a boundary crossing, excursion set model associated with the 
Poisson distribution in the initial Lagrangian space. The lin- 
ear barrier shape considered in this paper is a limiting form 
of that associated with the Poisson case. This correspon- 
dence is important, since the GPD is in good agreement 
with numerical simulations of clustering from Poisson ini- 
tial conditions (Itoh, Inagaki & Saslaw 1993 and references 
therein), provided the variance is allowed to be a free func- 
tion of scale. Bouchet & Hernquist (1992) showed that the 
GPD also describes clustering from white-noise initial con- 
ditions well. This just reflects the (intuitively obvious) fact 
that clustering from Poisson and white-noise initial condi- 
tions should evolve similarly. Since the Inverse Gaussian is 
just a limiting form of the GPD (just as the white-noise 
Gaussian is a limiting form of the Poisson distribution), the 
Inverse Gaussian should provide a good approximation to 
the Eulerian distribution measured in simulations of cluster- 
ing from white noise initial conditions, in the regime where 
discreteness effects are unimportant. Fig. ^ shows this ex- 
plicitly. 

The simulations of clustering from white noise initial 
conditions used here are the same as those studied by Mo 
& White (1996), where they are described in more detail. 
The simulations contain 10^ particles in a cube with peri- 
odic boundary conditions. The figure shows the evolved Eu- 
lerian distribution function for a few representative choices 
of the comoving scale. Each panel shows results for counts in 
spherical cells with comoving radii 0.02, 0.04, 0.08, and 0.16 
times the box size. For each cell size, the thin curves show 
the distribution of counts averaged over 30, 000 cells. The 
two panels show results at two different expansion factors a, 
where a = 1 initially. 

Comparison with the linear barrier model developed 
here is a little tricky. In principle, the model requires that 
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Figure 3. The Eulerian distribution function for a range of comoving scales, shown at two different output times, labelled by the 
expansion factor a (initially a = 1). Thin curves show the distribution measured in simulations of clustering from white noise Gaussian 
initial conditions. Thick curves show Inverse Gaussian distributions that have the same variance as the thin curves. Larger cells have 
narrower distribution functions. 



the free parameter of the Inverse Gaussian distribution 
(equation should scale as oc 1/V. In fact, although 
in the simulations is oc 1/V initially, it scales differently 
at later times (e.g. Hamilton et al. 1991; Jain, Mo & White 
1995). Therefore, when comparing equation (^o|) with the 
simulations, on any scale V, we assume that ^m(V) is given 
by the value measured in the simulations. Since is re- 
lated to the variance, the thick curves in Fig. ^ show In- 
verse Gaussian distributions that have the same variance as 
the distributions shown by the thin curves. When this is 
done, the Inverse Gaussian (equation describes the sim- 
ulation results reasonably well on all scales at all times. In 
fact, the Inverse Gaussian fits p{S) at least as well as the 
Press-Schechter mass function (equation fits the distri- 
bution of cluster masses (see, e.g., Lacey & Cole 1994). Since 
the evolution of the variance is reasonably well understood 
(Hamilton et al. 1991; Jain, Mo & White 1995) the results 
for the linear barrier presented here may be of more than 
just academic interest. 

For example, in the linear barrier model, when haloes 
first viri alize, they are uncorrelated with each other (Sec- 
tion 3.4). What correlations do exist arise primarily from 



volume exclusion effects. That is, in the initial Lagrangian 
space, haloes occupy a volume that is proportional to their 
mass. Since no other halo can occupy the region taken up 
by a given halo, this exclusion effect gives rise to (anti)- 
correlations, at least on small scales. This is probably a 
generic consequence of the white-noise initial conditions, and 
not specific to the linear barrier model. 



4 DISCUSSION AND EXTENSIONS 
4.1 Summary of the model 

This paper develops a model that allows one to provide an 
approximate description of the spatial distribution of dark 
matter, as well as dark matter haloes, in a self-consistent 
way. It can also be used to estimate the density profiles 
of haloes, and to quantify how matter flows in and out of 
Eulerian cells. The model is described in Section ^ 

The model assumes that, at the epoch z, the comov- 
ing Eulerian size R of an initial Lagrangian region Rq, is 
related to its initial overdensity So- This assumption is mo- 
tivated by the spherical collapse model. The additional as- 
sumption that the number of such regions is conserved, only 
their comoving size changes, allows one to compute statis- 
tics in the Eulerian space by taking appropriate averages 
over the relevant Lagrangian space quantities. Since the ini- 
tial Lagrangian distribution can be computed (e.g. Mo & 
White 1996; Sheth & Lemson 1998), so can the Eulerian 
distribution. The algorithm for taking these appropriate av- 
erages is related to the solution of an excursion set, barrier 
crossing problem, in which the barrier shape is given by the 
spherical collapse model. 

A constant barrier is a special case of the barrier stud- 
ied here. The first crossing distribution associated with a 
constant barrier can be used to estimate the distribution of 
virialized halo masses (Bond et al. 1991). In the context of 
this paper, this special case corresponds to the limit of van- 
ishing Eulerian cell size; the halo mass function is simply 
the Eulerian probability distribution function in the limit of 
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vanishing Eulerian cell size. Thus, the model here shows how 
the Bond et al. (1991) construction that yields the halo mass 
function can be extended to yield the Eulerian probability 
distribution function. The results of Mo & White (1996) are 
also easily understood within the context of this model. 

In Section ^ the linear barrier was used to show ex- 
plicitly how the model works. This was done for a number 
of reasons. Firstly, the spherical collapse barrier problem 
must be solved numerically, whereas, for the linear barrier, 
most interesting quantities can be computed analytically. 
Secondly, within the context of this model, the linear bar- 
rier can be understood as arising from a simple variant of 
the spherical collapse model (Section 3.1). For this barrier 
shape, the as soci ated Eulerian distribution is Inverse Gaus- 
sian (Section 3.3), and this is good agreement with numerical 
simulations of clustering from white noise initial conditions 
(Fig. This provides yet another reason for considering 
the linear barrier. 

In the model, correlations between haloes and mass 
arise from the simple fact that a cell which is known to 
contain a halo certainly contains that halo's mass. Correla- 
tions between haloes arise primarily from volume exclusion 
effects. That is, in the initial Lagrangian space, haloes oc- 
cupy a volume that is proportional to their mass. Since no 
other halo can occupy the region taken up by a given halo, 
this exclusion effect gives rise to (anti)-correlations, at least 
on small scales. This is probably a generic consequence of the 
white-noise initial conditions, and not specific to the linear 
barrier. 

For the linear barrier model, the evolution of halo-halo 



correlations can be written analytically (Section 3.4). Sec- 



tion 3.5 showed that the two-linear-barrier problem can also 



be solved analytically; in the model, this can be used to 
provide information about the density run around random 
positions in Eulerian space. Although the necessary formu- 
lae were all derived there, we did not pursue this further. 



4.2 Scaling properties of the model 

The Inverse Gaussian distribution was derived here by con- 
sidering the distribution of crossings of a linear barrier by 
Brownian motion. The relation between this barrier shape 
(equation ^ and that required by the spherical collapse 
model (equation was discussed in Section 3.1. 

The thick solid lines in Fig. H show this spherical 



collapse barrier (equation |l^ for representative values of 
V oc R^. Curves for larger R look just like those for smaller 
R except that are shifted to the left (since the a::-axis of 
the plot shows the negative of the logarithm of Vb). Curves 
for the same comoving R but at higher redshift are sim- 
ply multiplied by {1 + z), so the zero crossing remains the 
same. Dashed curves show the corresponding quantities for 
the linear barrier used in Section ^. This figure shows that 
changing the comoving size R at fixed z is equivalent to 
rescaling the Sq axis, and changing z at fixed comoving ra- 
dius R is equivalent to rescaling the So axis. This provides 
a strong constraint on the form of the distribution of first 
crossings. 

Namely, it implies that f{So,5o\R) should have a self- 
similar form, in that it should be a function of two parame- 
ters, one which is related to the scaling of the axis, and 
another which relates the scaling of the 5*0 axis to that of 




Figure 4. Dependence of the barrier shape on comoving Eule- 
rian size R and redshift z. Solid curves show B = 5o(Ro\R, z) of 
equation ([lo[), dashed curves show equation (po|). For white noise. 
So oc 1/Vo (X l/ijg. 



the So axis. Since the first crossing distribution is a function 
of only two parameters, the associated Eulerian distribution 
is also. It is convenient, then, to think of these parameters 
as being the mean and variance of the Eulerian distribution. 
This is an attractive feature of the model, since Colombi et 
al. (1997) find that, to a good approximation, the Eulerian 
distributions measured in TV-body simulations of clustering 
from white noise initial conditions are functions of the mean 
and the variance only. 



4.3 Scale— free initial conditions 

Initial Gaussian random fields with scale-free power spectra, 
P{k) oc fc" and —3 < n < 1, are a simp le g eneralization of 
white-noise (for which n — 0). Equation (hsf) shows that, for 
these spectra, 

SoocAfo"", where Q=(n + 3)/3, (55) 

which means that A = Mo/V = (Sv/So)^/°, where Sy oc 
V~" . This means that, when plotted as a function of {So, So), 
the shape of the spherical collapse barrier depends on the 
initial power spectrum. Although the barrier shape when 
n 7^ is different than for white noise, it still has the same 
scaling as the white noise barrier, since changing the Eule- 
rian scale R is still equivalent to a simple shift of the barrier 
along the log(5'o) axis. So, the barrier crossing and the as- 
sociated Eulerian distributions will still be functions of two 
parameters. This is in qualitative agreement with the results 
of Colombi et al. (1997). 

For initial conditions that are not scale-free, changing 
the Eulerian scale R is no longer equivalent to a simple shift 
of the barrier along the log(S'o) axis. Thus, the barrier cross- 
ing distribution is no longer self-similar, and the associated 
Eulerian distribution is no longer a function of just two pa- 
rameters. Of course, two parameter fits may remain a good 
approximation; at issue is the rate of change of the slope 
of the power spectrum over the scales at which most of the 
trajectories cross the barrier. 

While these are encouraging features of the model, ex- 
tending it to describe clustering from initial conditions that 



© 0000 RAS, MNRAS 000, 000-000 



An excursion set model 13 



are different from wliite noise is not straiglitforward. The 
problem is one of ensuring correct normalization. For scale- 
free initial conditions, equations (|5^) and (p^), and the re- 
quirement that the associated Eulerian distribution have 
unit mean (equation means that the first crossing dis- 
tributions must satisfy 



fiSo,5o\R)dSo = 



7^ 



f(So,So\R)dSo = 1. (56) 



The previous section showed that, for a white-noise Gaussian 
field, the distribution of first crossings of the linear barrier 
has an associated Eulerian distribution which has unit mean, 
as required. However, in general, not all barrier shapes are 
compatible with this requirement. 

For example, recall that equation ( |ll[ ) with Sco = 1.5 
is a good approximation to the spherical collapse model. 
When n = —1, a = 2/3, so the associated barrier, when 
written as a function of So, is the same (linear) as it was 
in Section ^ Now, for a linear barrier, the first passage dis- 
tribution is Inverse Gaussian. However, when a = 2/3 and 
f(So,5o\R) is Inverse Gaussian, then the second equality of 
equation (feq) is not satisfied. Thus, it appears likely that, 
for arbitrary power spectra, the excursion set approach de- 
veloped here will not work if at least some parametrizations 
of the spherical collapse barrier are used. 

This means that it may be more fruitful to consider the 
inverse problem to that considered in this paper. 



4.4 The inverse problem and the Generalized 
Inverse Gaussian distributions 

In this paper, the barrier shape was specified, then the first 
crossing distribution was obtained, and finally, the associ- 
ated Eulerian distribution was derived. It may be more use- 
ful to first specify an Eulerian distribution. Then, transform 
it to a first crossing distribution, from which the shape of 
the barrier can be inferred. This ensures that all normal- 
ization requirements have been satisfied. Once the barrier 
shape is known, the halo-mass and halo-halo correlations 
can be worked out just as they were in this paper. Finally, 
this barrier shape can be compared with that required by 
the spherical collapse model. 

In this context, it is worth mentioning that there is 
a class of distributions which provides convenient general- 
izations of the Inverse Gaussian, and may provide useful 
approximations to the Eulerian distributions measured in 
simulations of gravitational clustering from scale-free initial 
conditions. 

The Generalized Inverse Gaussian (GIG) distribution 
with index A, scale parameter rj, and concentration param- 
eter uj is 



fx{x\T],ijj)dx ■ 



2Kx{uj) 



da; 



a: > 0, 



(57) 



where K\ (u) is a modified Bessel function of the third kind 
with index A. For a given A, these are functions of just two 
parameters, 77 and oj, so they have the self-similar property 
discussed above. For these modified Bessel functions. 




K\[ulj) — K-\{uj) and 
Kx+i{lu) = 2{X/lu)Kx{lo) + Kx-i{lo). 



(58) 



Variance 

Figure 5. 83 and 5*4 versus variance for the Generalized Inverse 
Gaussian distributions associated with initial power spectra with 
slopes n = (filled circles), —1 (open circles), and —2 (stars). The 
soHd curves show the values measured in numerical simulations 
of clustering from these initial conditions, as parametrized by 
Colombi et al. (1997). Dashed curves show their estimates of the 
allowed range of values. 



Therefore, when A = — (2q)^^, then the GIG distributions, 
with X = So, rj = Sv and lu = S^o/^v satisfy equation (^jj), 
so the associated Eulerian distributions are correctly nor- 
malized and have unit mean. When n = 0, then a = 1, so 
A = —1/2, and this distribution is the same as the Inverse 
Gaussian distribution considered in section ^. 

The symbols in Fig. |^ show the values of S3 and 54 
(defined similarly to equation ^) plotted versus variance for 
these GIG distributions when n = (filled circles), —1 (open 
circles), and —2 (stars). The curves were computed using the 
extended perturbation theory fitting functions of Colombi et 
al. (1997). The solid curves show their fits to the values of S3 
and S4 measured in numerical simulations of clustering from 
the corresponding scale-free initial conditions, and dashed 
curves give an estimate of the allowed range of values. While 
the GIG values do not fit the numerical results, they do show 
similar trends. 

The problem, then, is to determine the shape of the 
barrier for which the GIG distribution is the first crossing 
distribution, and to then compute the associated halo dis- 
tribution. The solution of this inverse problem is the sub- 
ject of ongoing work, where the barriers associated with the 
distributions given by the extended perturbation theory of 
Colombi et al. (1997) are also used as the Eulerian distribu- 
tions. 
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APPENDIX A: THE FIRST PASSAGE TIME 
AND THE LINEAR BARRIER 

This Appendix presents a derivation of the first time S that a 
particle undergoing Brownian motion with zero drift reaches 
the linear barrier 



(Al) 



having started from the origin {S,S) — (0,0). Let fu{S) de- 
note the first passage time distribution associated with this 
barrier. The subscript v indicates the height of the barrier 
at 5" = 0. 

The assumption of zero drift means that, in the absence 
of the absorbing barrier, the mean distance of a particle from 
the S axis, averaged over the ensemble of Brownian walks of 
the type shown in Fig. |l| is zero. To be more precise, recall 
that equation ( [27| ) gives the probability p{S, So) d6o that a 
trajectory has value between So and So + dSo, at 5*. So, for 
those trajectories, 

(So) = / Sop{S,So) dSo = for aU S, (A2) 



and 



(So) = y 5o pis, So) dSo = S + (So)^ = S at S. (A3) 

Since (So) ~ for all 5", these are said to be Brownian walks 
with zero drift. 

Consider another linear barrier which is a simple shift 
of fj, along the S axis of the barrier above: 



Suj{S) = ^ + So{S) = ^ -\- V — 13 S = iu — f3 S. 



(A4) 



Let fui{S) denote the first passage time distribution associ- 
ated with this barrier. The subscript uj indicates the height 
of the barrier when 5* = 0. 
Then 



f^is)ds= / ds' Ms'),u{s\s')ds, 



(A5) 



where fuii^{S\S') denotes the first passage time distribution 
to Sij{S), given that the particle started at [S' , S,^{S')] in- 
stead of the origin. Since the particle is undergoing Brow- 
nian motion, fui,y{S\S') is the same as for a particle that 
starts at the origin, but sees a barrier 

S^{S) = S^{S) - S,{S') =yi~l3{S^ S'). (A6) 

So, we can write 

/^.(5|S") =/m(5-S"), (A7) 

where the subscript /i indicates that the height of the barrier 
above the starting position at S' is fj,. 

This means that equation ( |A5| ) is a convolution equa- 
tion: 

US) = f^+,.{S) = / US') us - S') dS'. (A8) 
Jo 

When /3 = 0, then the barrier is constant, and the first 
passage distribution is known to be 

1/2 



US)dS = 



( — 

\2nS 



-i/^/2S 



dS 
S ' 



(A9) 



It is straightforward to verify that this expression satisfies 
the convolution equation above. Therefore, the solution to 
the convolution equation when /3 7^ must be 



US)dS- 



\2nS 



1/2 



-{iy-l3Sy/2S 



dS 

s • 



(AlO) 



A more formal derivation can be found in, e.g., Schrodinger 
(1915) or Cox & Miller (1967). The following derivation fol- 
lows that in Kao (1996) closely. 



First, multiply both sides of equation ( A5) by exp(— tS), 
and integrate over all S. This gives the Laplace transform 
C{uj,t) of fu}{S). Equation ([AE|) implies that 



C(uj,t) = I dS e 

Jo 



/ ds' us')Uis\s') 

dS' US') I dS US-S') e"*^ 



dS' US')e''"^' / dS" US")e-''^" 
Jo 
= C{u,t) C{fi,t), (All) 

which is the same as 

C{u + ii.t) ^ C{v,t) C{ii,i). (A12) 

This is a functional equation with solution Ciy) = , 
where C > is some constant, and we have not bothered to 
write the dummy variable t. The problem, then, is to find 
C. 

To do so, notice that the crossing of a barrier which de- 
creases as S increases, by Brownian walks with zero drift, is 
equivalent to the crossing of a constant barrier by trajecto- 
ries that, in the mean, drift upwards as S increases, provided 
that the drift is chosen correctly. This correct choice of drift 
corresponds to choosing 

{So)=l3S at 5". (A13) 
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Of course, the variance remains the same, so that 



at S. 



(AM) 



So the problem of the crossing of the hnear barrier = 
y — pS by walks with zero drift is transformed to the problem 
of the crossing of a constant barrier of heig ht v b y walks with 
upward drift (5. In particular, the relation {^.12) remains the 
same. 

So, to find C, assume that the probability of first cross- 
ing the barrier u in the first small increment s along the 
S axis is negligible. Suppose that at S = s, the random 
trajectory has value 5o < v. Then 



MS) = 

so that 



f^{S\So,s)p{s,So) dSo 
f^-So(S - s)p(s,So) dSo, 



(A15) 



dS e-*^ J f^^so {S - s) p(s, So) dSo 
dSop{s,5o) / dS U-s„iS~s)e-"^ 

■J s 

dSop{s,5o) C{u-So,t). (A16) 



In what follows, we will sometimes omit the dummy variable 
t: C{u) = jC{i',t). As s and So are small. 



Ciu - So) ~ £M - So + ^C"{i^), 
and 

These expressions imply that 



(A17) 
(A18) 



^ (^i^ts) (^ciu)-{5o)c'{u)+^-^/:"{u) 

= {l-ts)(^C{u)-(3sC'{u)+'-±^L"{u)^ 
^ C{i^)-tsC{u)-/3s£'{u) + ^C"{u), (A19) 



The second line follows from (A13) and (A14), and the third 
is correct to lowest order in s. Dividing by s, and taking the 
limit s —> 0, reduces this to 



t£M = -/3£'M + -£"H. 
However, 

= e"*^"^, so 
These expressions in ( |A2C| ) imply that 

+ 2PC - 2t = 0, 
so that, if C > 0, then 
C = ~f3+ + 2t, 
since /3 and t axe both positive. Thus, 



(A20) 



Ciu,t) 



Ce ^'^ , and C" (u) = C^e 



(A21) 
(A22) 
(A23) 



Inverting this Laplace transform gives the Inverse Gaussian 
distribution. 

This problem can also be formulated in the context 
of the diffusion equation. In the notation of Lacey & Cole 
(1993), the associated diffusion equation is 



dS ^ ' dS 2dS^' 



(A24) 



where /3 is the drift term, and Q{S, S, v) represents the solu- 
tion to this equation in the presence of an absorbing bound- 
ary &t S = u. Bond et al. (1991) and Lacey & Cole (1993) 
considered the case /3 = 0. 
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